Automated method for identifying landmarks within an image of the brain

ABSTRACT

A method is disclosed for obtaining the location of a landmark in an MR image of a brain. In a first step, a region of interest in a plane within the MR image containing the landmark is defined. In a second step, the ROI is binarised into foreground and background voxels based on at least one threshold selected using anatomical knowledge. In a third step a set of object voxels is identified from the foreground voxels, excluding voxels which were only classified as object due to proximity of cortical and non-cortical structures. This can be done by morphological processing which reclassifies voxels which may have been incorrectly classified as object, followed by restoring voxels due to the partial volume effect and/or morphological erosion/opening. In a fourth step, an automatic process is then carried out to identify one or more landmarks in the modified binarised image.

FIELD OF THE INVENTION

The present invention relates to methods for automatically identifyinglandmarks within images of a brain.

BACKGROUND OF INVENTION

The Talairach transformation is widely used for analysis of neurologicalimages. It involves identifying eight landmarks, which are used todefine a coordinate system. The Tailairach landmarks subdivide the braininto 12 cuboids, and the Tailairach transformation is to warp the imageswithin each cuboid linearly. In this way the brain images are normalisedby a three-dimensional piece-wise linear warping. This scheme hasseveral applications, in particular because it makes it possible tocompare neurological images from different individuals. One improvementon this scheme, while following its conceptual rationale, is theimprovement of the definitions of the landmarks, to give “modifiedTalairach landmarks” (as defined in the article “Modified TalairachLandmarks”, W. L. Nowinski, Acta Neurochirurgica, 2001, 143, p1045-1057, the disclosure of which is incorporated herein by reference).

A principal advantage of the Talairach transform is its simplicity.Although numerous non-linear image registration techniques are known, inprinciple providing a higher accuracy, the non-linear techniques havelimitations which make them difficult to use beyond a researchenvironment. In particular, a prohibitively high computational time isrequired. Whereas the Talairach transformation can be performed in lessthan a second in a standard personal computer, some non-linear methodsrequire days of computation. More fundamentally, the non-linear methodsare conceptually complex, and must be treated as “black boxes”, whichlimits their clinical acceptance.

Furthermore, there is no established methodology for validation of thenon-linear registration techniques.

One drawback of the Talairach transformation, however, is that theidentification of the landmarks has not so far been automated reliably,so that when using conventional software which employs the Talairachtransformation the landmarks still have to be identified by userinteraction (e.g. R. W. Cox, “AFNI: Software for Analsysis andVisualization of Functional Magnetic Resonance Neuroimages”, Computerand Biomedical Research, 1996, 29, p 162-173). Quite apart from the timethe interactive identification of the landmarks takes, differentindividuals are liable to locate the landmarks in slightly differentpositions, which reduces the robustness of the method.

SUMMARY OF THE INVENTION

The present invention aims to provide new and useful techniques andapparatus for identifying automatically landmarks within neurologicalimages, and in particular the Talairach landmarks or modified Talairachlandmarks.

The invention makes use of the concepts “foreground voxels” and“background voxels”. “Binarising” a set of voxels means to divide thevoxels into these two categories.

Furthermore, certain voxels are “object voxels”. Object voxels arevoxels in a location which corresponds to a physical identity (cerebralstructure), such as, in the context of this invention, the union of greymatter and white matter.

In general terms, the invention proposes that the location of one ormore landmarks is obtained by the following steps:

(a) a region of interest in a plane within an MR image of a brain andcontaining the landmark(s), is identified;

(b) the region of interest is binarised into foreground and backgroundvoxels based on at least one threshold selected using anatomicalknowledge,

(c) a set of object voxels is identified from the foreground voxels,excluding voxels which were only classified as foreground voxels due toproximity of cortical and non-cortical structures,

(d) object voxels are identified from the background voxels due to thepartial volume effect and morphological erosion/opening, and

(e) an automatic process is then carried out to identify one or morelandmarks in the modified binarised image.

The step of reclassifying voxels may be performed in two stages: (i) afirst stage of morphological processing, and (ii) a second step ofrestoring voxels which have been incorrectly reclassified in themorphological processing.

Anatomical knowledge may also be used in the reclassifying step, e.g. byusing the expected shapes of cortical and/or non-cortical structures tomodify the classification.

The threshold is preferably selected by the steps of:

(i) using prior knowledge about the image to derive an intensity rangeof voxels in the said region of interest;

(ii) obtaining a frequency distribution of the intensities within thesaid intensity range of voxels within the said region of interest; and

(iii) using the said frequency distribution to derive the threshold.

The intensity threshold may be selected by minimising a function whichis a sum of the variances of the intensities below and above thethreshold.

Optionally, this function may be a weighted sum defined based on twoconstants W₁ and W₂. This is referred to here as a “range constrainedweighted variance method”.

For example, labelling the possible values of voxel intensity by aninteger index i and their respective frequencies by h(i), and writingthe lower and upper intensities respectively as r_(low) and r_(high),the weighted sum is given by${{\theta_{RCLWV}\left( {W_{1},W_{2}} \right)} = {\max\limits_{r_{k}}\left( {{{\Pr\left( C_{1} \right)}{D\left( C_{1} \right)}W_{1}} + {{\Pr\left( C_{2} \right)}{D\left( C_{2} \right)}W_{2}}} \right)}},$where Pr(·) denotes the class probability$\left( {{\Pr\left( C_{1} \right)} = {{\sum\limits_{i = r_{low}}^{r_{k}}\quad{{h(i)}\quad{and}\quad{\Pr\left( C_{2} \right)}}} = {\sum\limits_{i = {r_{k} + 1}}^{r_{high}}\quad{h(i)}}}} \right),,$and D(C₁) and D(C₂) are given by:D(C ₁)=(μ₀−μ_(T))² and D(C ₂)=(μ₁−μ_(T))², where${\mu_{T} = {\sum\limits_{i = r_{low}}^{r_{high}}\quad{i \times {h(i)}}}},{\mu_{0} = {\sum\limits_{i = r_{low}}^{r_{k}}\quad{i \times {h(i)}}}}$and$\mu_{1} = {\sum\limits_{i = {r_{k} + 1}}^{r_{high}}\quad{i \times {{h(i)}.}}}$

The steps (a) to (d) may be performed repeatedly for differentlandmarks. One or more landmarks may be located for each of the startingplanes of the MR image.

In the case of identifying the I landmark (which is a position on acoronal plane) the left and right halves of the brain may be treatedpartly separately, such that the I landmark is identified using the halfof the brain which appears to have been segmented more accurately.

The invention may be expressed as a computer-implemented method, oralternatively as a computer system arranged to perform such a method.

BRIEF DESCRIPTION OF THE FIGURES

Preferred features of the invention will now be described, for the sakeof illustration only, with reference to the following figures in which:

FIG. 1 is a flow chart of a method which is an embodiment of theinvention;

FIG. 2 is a flow chart showing in more detail sub-steps of step 3 of theflow chart of FIG. 1;

FIG. 3, which is composed of FIGS. 3(a) to 3(d), shows an example ofperforming the steps to identify the A, P, L and R landmarks in themethod of FIG. 1;

FIG. 4, which is composed of FIGS. 4(a) to 4(d), shows an example ofperforming the steps to identify the S landmark in the method of FIG. 1;and

FIG. 5, which is composed of FIGS. 5(a) to 5(d), shows an example ofperforming the steps to identify the I landmark in the method of FIG. 1

DETAILED DESCRIPTION OF THE EMBODIMENTS

Referring firstly to FIG. 1, the steps of a method which is anembodiment of the invention are shown.

In step 1, a dataset which is a neuroimage (i.e. an image of a brain) isloaded. This image is a three dimensional data, typically obtained froma number of scans in different respective planes.

From this data, the midsagittal plane (MSP) is determined. This ispreferably done using the method disclosed in WO02/069827, “Method andapparatus for determining symmetry in 2D and 3D images”, by Hu andNowinski (the disclosure of which is incorporated herein by reference).However, the invention is not limited in this respect, and any othertechnique for determining the MSP may also be applied. Indeed, it wouldalso be possible within the scope of the invention for the input data tospecify the MSP.

The coordinates of the anterior commissure (AC) and posterior commissure(PC) are then determined automatically. This can be done by the methoddisclosed in WO02/43003, “Methods and apparatus for processing medicalimages”, by Nowinski and Thirunavuukarasuu (the disclosure of which isincorporated herein by reference), although once more the invention isnot limited in this respect.

In step 2 of FIG. 1, the data is normalised to occupy a predefinedgray-scale range. According to the standard radiological convention, wewrite the coordinate system as (x,y,z), where the x-axis is from thesubject's right to left, the y-axis is from anterior to posterior, and zfrom superior to inferior. Thus, an xz-plane is a coronal plane, ayz-plane is a sagittal plane, and an xy-plane is an axial plane. Letg(x,y,z) denote the gray level of the input data at a voxel at position(x,y,z), and let g₀ and g₁ denote the grey levels such that one percentof the voxels have an intensity less than g₀ and one percent of thevoxels have an intensity greater than g₁. Then, we obtain a normalisedgray level ĝ(x,y,z)which, for a given position (x,y,z) is given by:$\begin{matrix}{{\hat{g}\left( {x,y,z} \right)} = 0} & {{{if}\quad{g\left( {x,y,z} \right)}} \leq g_{0}} \\{{\hat{g}\left( {x,y,z} \right)} = \frac{{g\left( {x,y,z} \right)} - g_{0}}{g_{1} - g_{0}}} & {{{if}\quad g_{0}} < {g\left( {x,y,z} \right)} < g_{1}} \\{{\hat{g}\left( {x,y,z} \right)} = 255} & {{{if}\quad{g\left( {x,y,z} \right)}} \geq g_{1}}\end{matrix}$

Each slice of the normalised data has its own co-ordinate system (u,v)where u is the horizontal direction and v is the vertical direction.

In step 3, the position of the A, P, L and R landmarks is located. Thisis done by the series of steps shown in FIG. 2.

Firstly, in step 3.1 a region of interest (ROI) is defined. These may bethe voxels within the skull. These voxels can be determined by thefollowing steps:

(a) A histogram-based thresholding method is used to binarise the APplane (as used for example in M. E. Brummer, R. M. Mersereau, R. L.Eisner, R. J. Lewine, “Automatic detection of brain contours in MRI datasets”, IEEE Transactions on Medical Imaging 1993; 12(2), p 153-166).

(b) A morphological closing operation is performed using a 3×3structuring element (SE) to connect small gaps in the ROI.

(c) The largest connected component is identified, and the holes withinthe component are filled.

FIG. 3(a) shows the AP plane in a typical example of the use of thismethod. FIG. 3(b) shows the corresponding ROI.

In step 3.2 an optimum threshold is determined, based on therange-constrained weighted variance thresholding method. This includesfollowing steps, which are explained in a separate patent application bytwo of the present inventors: “Methods and apparatus for binarizingimages”, Singapore patent application 200307531-4, by Q. M. Hu, Z. Hou,and W. L. Nowinski, which was still unpublished at the priority data ofthe present application, and the disclosure of which is incorporated byreference.

Firstly, prior knowledge of the image is used to define an ROI which isa subset of the image. This process can be done by whatever means,either automatic, semi-automatic, or even manual. Then an analysis isperformed on the frequency of occurrence of intensities within the ROI,and a range of frequencies is defined, again using prior knowledge.

For example, without losing generality, we denote the image to beprocessed as f(x), where f(x) is the gray level at a voxel labelled x.It is further supposed that the processed image has L gray levelsdenoted by r_(i) where i is an integer in the range 0 to L-1 and r₀<r₁<. . . r_(L-1). It is also assumed that the object of interest has higherintensity values than the background. Suppose that due to priorknowledge or test we know that the proportion of the ROI which isoccupied by the object is in the percentage range per₀ to per₁.

Let h(i) denote the frequency of gray level r_(i), and let H(i) denotethe cumulative frequency which is${\sum\limits_{i^{\prime} = 0}^{i}\quad{h\left( i^{\prime} \right)}},$where i′ is an integer dummy index. Considering two values of i writtenas m and n, the frequency of intensities in the range r_(m) to r_(n) is$\sum\limits_{i^{\prime} = m}^{n}\quad{{h\left( i^{\prime} \right)}.}$Thus, we can use peroto calculate a gray level r_(low), such that we aresure that all the voxels having lower intensity represent background.r_(low) can be written as: $\begin{matrix}{r_{low} = {\min\limits_{i}{\left\{ i \middle| {{H(i)} \geq {per}_{0}} \right\}.}}} & (1)\end{matrix}$

Similarly, we can use per₁ to calculate a gray level r_(high) such thatwe are sure that all voxels having higher intensity represent theobject: $\begin{matrix}{r_{high} = {\min\limits_{i}{\left\{ i \middle| {{H(i)} \geq {per}_{1}} \right\}.}}} & (2)\end{matrix}$

Let r_(k) fall within the range r_(low) to r_(high), and suppose thatthe voxels of the ROI are in two classes C₁ and C₂, where C₁ is voxelsof the background class and consists of voxels with gray levels r_(low)to r_(k), and C₂ is voxels of the object class and is composed of voxelswith gray levels r_(k+1) to r_(high). The range-constrained weightedvariance thresholding method maximises the “weighted between-classvariance” defined as:${{\theta_{RCWV}\left( {W_{1},W_{2}} \right)} = {\max\limits_{r_{k}}\left( {{{\Pr\left( C_{1} \right)}{D\left( C_{1} \right)}W_{1}} + {{\Pr\left( C_{2} \right)}{D\left( C_{2} \right)}W_{2}}} \right)}},$where W₁ and W₂ are two positive constants selected by the user andrepresenting the weights of the two respective class variances, Pr(·)denotes the class probability, i.e.${{\Pr\left( C_{1} \right)} = {\sum\limits_{i = r_{low}}^{r_{k}}\quad{h(i)}}},{{\Pr\left( C_{2} \right)} = {\sum\limits_{i = {r_{k} + 1}}^{r_{high}}\quad{h(i)}}},$and D(C₁) and D(C₂) are given by:D(C ₁)=(μ₀ −μ _(T))² and D(C ₂)=(μ₁ −μ _(T))², where${\mu_{T} = {\sum\limits_{i = r_{low}}^{r_{high}}{i \times {h(i)}}}},{\mu_{0} = {{\sum\limits_{i = r_{low}}^{r_{k}}{{i \times {h(i)}}\quad{and}\quad\mu_{1}}} = {\sum\limits_{i = {r_{k} + 1}}^{r_{high}}{{i \times {h(i)}}.}}}}$

When W₁ is bigger than W₂, background homogeneity is emphasised.

Step 3.2 may be done by specifying per₀ and per₁ to be 14% and 28%respectively. The optimum threshold is denoted as θ₁.

In step 3.3, we segment the AP plane by assigning voxels to foregroundif they are bigger than θ₁, and otherwise assigning them to background.The binarised image is denoted as BWAP1 (u,v).

In steps 3.4 we reclassify the voxels: firstly by a morphologicalprocessing, then processing using anatomical knowledge, and finallyperforming a restoring step.

The sub-steps of the morphological processing and processing usinganatomical knowledge are as follows:

(a) A distance transform of the ROI is performed using the 2-3 metric(in this metric the distances between any two voxels is determined bydefining the shortest path of voxels between them, and adding thedistance increments along this path. The distance increments between twovoxels which are nearest neighbours in a direction parallel to one ofthe axes (4-connected nearest neighbours) is taken as 2, and thedistance between two voxels which are nearest neighbours diagonally(8-connected nearest neighbours) is taken as 3) and the distance codesare converted into distance indices by the method of Hu Q. M., “Two andthree-dimensional skeletonization”, WO 02/058008). The maximum distanceindex is denoted as maxDSkull.

(b) A morphological opening operation is performed with a 3×3 SE withrespect to BWAP1 (u, v), to obtain BWAP2(u,v).

(c) A morphological opening operation is performed with a 5×5 SE withrespect to BWAP1(u,v) to obtain BWAP3(u,v).

(d) The foreground components of BWAP3(u,v) are found. For eachforeground component, its minimum and maximum distance indices aredenoted as minD and maxD respectively. A foreground component is treatedas an object component when maxD-minD is bigger than a value (e.g. 20)or maxD is bigger than a second value (e.g. maxDSkull/2).

(e) The object voxels are excluded from the foreground voxels ofBWAP2(u,v). The connected foreground components of BWAP2(u,v) are found.The number of voxels of each foreground component are denoted asnosVoxel. A foreground component of BWAP2(u,v) is categorised as anobject component only when the shape of the component is not similar tomeninges. According to anatomical knowledge, meninges have a shapesimilar to the skull and are quite thin. So, when maxD-minD is smallerthan 0.1*nosVoxel, the component is highly likely to be a meninges.Otherwise, it is classified as an object component.

The lost object voxels are restored by the following steps:

(a) Object voxels far from the skull lost due to the morphologicalopening operation are restored. This is achieved by checking thenon-object voxels with a distance index greater than maxDSkull/4. Iftheir gray level is bigger than θ₁, the voxels are reclassified asobject voxels.

(b) Object voxels around the boundaries lost due to the morphologicalopening operation are restored. For each object boundary voxel (anobject boundary voxel is an object voxel having at least one non-objectvoxel as one of its 8 immediate neighbours), each of its 8 immediateneighbours is reclassified as an object voxel if its grey level isgreater than θ₁.

Note that this restoration is not performed around the most anterior(i.e. minimum v) and most posterior (i.e. maximum v) parts of thestraight line connecting AC and PC, to avoid the sinus and meninges.Specifically, suppose the maximum and minimum v coordinates of objectvoxels are minVap and maxVap respectively, and denote the coordinates ofAC on the AP plane as (acUap,acVap) and the coordinates of the PC on theAP plane as (pcUap, pcVap). The restoration is not carried out in thefollowing two regions:|u-acUap|<10 mm and |v-minVap|<3 mm.   (3)|u-pcVap|<10 mm and |v-maxVap|<3 mm,   (4)where |x| stands for the absolute value of x.

(c) Object voxels (that is, all voxels which are physically either GM orWM) lost due to the partial volume effect are restored. Since thestatistics of GM, WM, CSF, air, meninges and bones are not available,the partial volume effect is alleviated by reducing the threshold by 10.That is, for any object boundary voxel, each of its immediate 8neighbours is re-classified as an object voxel if its gray level isbigger than (θ₁−10). The restoration is not carried out in the tworegions defined by formulae (3) and (4).

In step 3.5, the minimum and maximum v coordinates of the object voxelsare taken respectively as the v coordinates of the A and P landmarksrespectively. Similarly, the minimum and maximum u coordinates of theobject voxels are taken as the u coordinates of the L and R landmarksrespectively. Note that the u-coordinate in the AP plane corresponds tothe x-coordinate of the three-dimensional volume, and the v-coordinatein the AP plane corresponds to the y-coordinate in the three-dimensionalvolume.

FIGS. 3(c) and 3(d) show the segmented AP plane, and the 4 landmarksoverlaid over it. The white horizontal lines show the v coordinates ofthe A and P landmarks, while the white vertical lines shown the ucoordinates of the L and R landmarks on the AP plane.

In step 4 of FIG. 1, we determine the position of the S landmark. TheVPC plane is a coronal slice perpendicular to both the MSP and the APplane, and it passes through the PC. To determine the position of the Slandmark, we only need to determine its v co-ordinate. The v coordinateof the S landmark is the smallest v coordinate of all the corticalvoxels on the VPC plane. The S landmark is localized by segmentation ofa virtual slice aVPC(u,v) with the close skull.

The VPC plane is denoted as VPC(u,v), the coordinates of the PC withinVPC(u,v) are denoted as (pcU, pcV).

In step 4.1, a virtual plane aVPC(u,v) is constructed in the followingway:

(a) aVPC(u,v)=VPC(u,v) when v is not bigger than pcV.

(b) aVPC(u,v)=VPC(u,2pcV-v) when v is bigger than pcV and smaller than2pcV.

FIGS. 4(a) and 4(b) show a VPC and the corresponding virtual slice aVPC.

In step 4.2, the ROI corresponding to aVPC is found. This procedure isdone in the same way in which the ROI for the AP plane was found above(step 4.1 above).

In step 4.3, the optimum threshold θ₂ is determined by therange-constrained weighted thresholding method, by specifying per₀ andper₁ to be 20% and 40% respectively.

In step 4.4, the aVPC plane is segmented using the optimum threshold θ₂,by the same set of sub-steps explained in step 4.3 above.

In step 4.5, the threshold θ₂ is adjusted using anatomical knowledge.Since the vertical line passing through the PC in the vicinity of the Slandmark should be interhemispheric fissure voxels, θ₂ should be higherthan the gray levels of voxels on the vertical line segment startingfrom 2 mm above and ending 2 mm below the object voxels with the minimumv coordinate found in step 4.4. If θ₂ is indeed bigger than this value,then it is not changed. However, if it is lower, it is modified to be 5plus the maximum gray level of the line segment, and the aVPC(u,v) isre-segmented with the modified threshold θ₂, by the same sub-steps asthose used to segment the AP plane in step 4.3 above, including the samemorpological opening operations.

In step 4.6, lost object voxels are restored. This is done by thefollowing steps:

(a) Object voxels around the object boundaries lost due to themorphological opening operation are restored. This is done by,reclassifying an non-object voxel which has an object boundary voxel asa nearest neighbour and which has a gray level greater than θ₂.

Note, however, that to avoid the sinus/meninges, the restoration is notcarried out in the region defined by:|u-pcU|<10 mm and |v-minVavp|<3 mm   (5)where minVavp denotes the minimum v coordinate of object voxels.

(b) Object voxels lost due to the partial volume effect are restored.This is done by reclassifying an non-object voxel having an objectboundary voxel as a nearest neighbour and which has a gray level higherthan (θ₂−10). The restoration is not, however, carried out in the regiondefined by equation (5).

The segmented aVPC after restoration sub-steps (a) and (b) is shown inFIG. 4(c).

In step 4.7, the v coordinate of the S landmark is the mimimum vcoordinate of all object voxels in aVPC. The v coordinate of S in theaVPC plane is equal to its z coordinate in the full three-dimensionalvolume.

FIG. 4(c) shows the segmented aVPC slice, and FIG. 4(d) shows theoriginal VPC overlaid by a horizontal line indicating the z coordinate(or equivalently v coordinate) of the S landmark.

In step 5, the position of the I landmark is determined. The VAC planeis a coronal slice parallel to the VPC plane, and passes through the AC.Only the v coordinate of the I landmark is need, and it is the maximum vcoordinate of all the cortical voxels on the VAC plane. It is determinedby segmenting the VAC plane directly constrained by anatomicalknowledge. It is assumed that the maximum difference in the z-coordinatebetween the AC and the I landmark is no more than 50 mm. We denote thatcoordinates of the AC in the VAC(u,v) by (acU, acV).

Specifically, the I landmark is obtained by the following sub-steps:

In step 5.1, VAC(u,v) is binarised by assigning all voxels with graylevels bigger than θ₂ to be foreground voxels, and the rest asbackground voxels. The binarised image is denoted as BWVAC1(u,v).

In step 5.2, the region around the AC is connected to make subsequent“seeding” feasible. Here “seeding” means to find connected componentsfrom a specified voxel (“seed”) with all the voxels in the componentbeing of the same type (i.e. background, foreground or object). Theregion around the AC is connected by setting the BWVAC1(u,v) toforeground when |v-acV| is smaller than 3 mm.

Step 5.3 employs a vertical line passing through the AC to separate theVAC into left and right halves. The voxels on the vertical line with av-coordinate greater than acV+3 mm are set to background. That is,BWVAC1(u,v) is set to background when v is bigger than (acVac+3 mm) andu is equal to acU. This has the effect that there is foregroundseparation in the neck region.

In step 5.4, a morphological opening operation is performed onBWVAC1(u,v) using a 3×3 SE, to give BWVAC2(u,v). This operation breaksweak connections between the cortex and the skull and between the cortexand the neck.

In step 5.5, a morphological erosion operation is performed using a 3×3SE, to give BWVAC3(u,v). This operation further breaks the connectionsbetween the cortex and the skull and between the cortex and the neck.

In step 5.6, we seed from (acU, acV), to obtain the foregroundcomponent. Then, we perform a morphological dilation on the seededforeground component with a 3×3 SE, to obtain BWVAC4(u,v). The serialoperations (erosion followed by seeding and dilation) are intended tobreak strong connections between the cortex and the skull, and betweenthe cortex and the neck, while retaining the original shape of thecortex.

In step 5.7, the maximum value of v for which a foreground BWVAC4(u,v)voxel exists having u smaller than acU, is found, and denoted as maxVL.

In step 5.8, the maximum value of v for which a foreground BWVAC4(u,v)voxel exists having u not less than acU, is found, and denoted as maxVR.

In step 5.9, the left half of BWVAC4 (u,v) (i.e. the values of u smallerthan acU) is restored in two substeps if (maxVL-acV) is smaller than 50mm:

(a) Firstly, the effects of the morphological opening operation arecounteracted, by changing any background voxel which has at least aforeground voxel of BWVAC4(u,v) as one of its 8 immediate neighbours andwhich has a gray level in VAC(u,v) greater than θ₂ to foreground.

(b) Secondly, the effects of the partial volume effect are counteractedby restoring any background voxel which has a foreground boundary voxel(a foreground boundary voxel is a foreground voxel, within its 8immediate neighbours there is at least a background voxel) ofBWVAC4(u,v) as one of its 8 immediate neighbours and which has a graylevel in VAC(u,v) bigger than (θ₂−10) to foreground.

Similarly, the right half of BWVAC4(u,v) (i.e. the half with u not lessthan acU) is restored by two corresponding steps when maxVR-acV issmaller than 50 mm.

Object voxels are all foreground voxels in this case.

In step 5.10, if both (maxVL-acV) and (maxVR-acV) are smaller than 50mm, the v coordinate of the I landmark is obtained as the v component ofthe object voxel in BWVAC4(u,v) with the biggest value of v. If one of(maxVL-acV) or (maxVR-acV) is smaller than 50 mm and the other is not,the v-coordinate of the I landmark is obtained as the maximum vcoordinate of all the object voxels from the side (i.e. the left side orright side) for which the maximum v coordinate of the object voxels issmaller than (50+acV). If both of (maxVL-acV) and (maxVR-acV) are biggerthan 50 mm, the v coordinate of the I landmark Is the maximum vcoordinate of all the object voxels on the side (i.e. left or right)whose maximum object v coordinate is smaller. Note that the v coordinateof BWVCA(u,v) corresponds to the z coordinate of the dataset.

FIG. 5(a) shows the original VAC, FIG. 5(b) shows the binarised VAC(i.e. BWVAC1), FIG. 5(c) shows the processed foreground (BWVAC4), andFIG. 5(d) shows the v coordinate (or equivalently z coordinate) of the Ilandmark overlaid on the original VAC plane.

Note that the set of sub-steps performed in steps 4 and 5 can beconsidered as corresponding to those shown in FIG. 2 for finding the A,P, L and R landmarks. That is, a ROI is identified; a threshold isselected (or, in the case of step 5, the threshold used is that same onederived in step 4); a segmentation is performed; then a reclassificationis performed; and finally the landmarks are identified.

Although only a single embodiment of the invention has been described indetail, many variations are possible within the scope of the inventionas will be clear to a skilled reader.

1. A computer-based method for locating one or more landmarks using anMR image of a brain, the method including the following automatic steps:(a) identifying a region of interest (ROI) with a plane of the MR image,the plane containing the landmark(s); (b) binarising the plane of the MRimage into foreground and background voxels based on at least onethreshold selected using anatomical knowledge; (c) identifying a set ofobject voxels from the foreground voxels, the set of object voxelsexcluding voxels which were only classified as foreground voxels due toproximity of cortical and non-cortical structures; (d) identifyingobject voxels from the background voxels due to the partial volumeeffect and/or morphological erosion/opening; and (e) identifying the oneor more landmarks using the object voxel.
 2. A method according to claim1 in which the step of identifying the object voxels is performed in twostages: (i) morphological processing which excludes foreground voxelswhich may not be object voxels, and (ii) restoring voxels which havebeen incorrectly excluded in the morphological processing.
 3. A methodaccording to claim 2 in which the step of identifying the object voxelsfurther includes applying anatomical knowledge to identify the objectvoxels.
 4. A method according to claim 3 in which the anatomicalknowledge is knowledge about the expected shapes of cortical and/ornon-cortical structures.
 5. A method according to claim 1, wherein thethreshold is selected by the steps of: (i) using prior knowledge aboutthe image to derive an intensity range of voxels said region ofinterest; (ii) obtaining a frequency distribution of intensities withinthe said intensity range of voxels within said region of interest; and(iii) using the frequency distribution to derive an intensity threshold.6. A method according to claim 5 in which the intensity threshold isselected by minimizing a function which is a sum of variances of theintensities below and above the threshold.
 7. A method according toclaim 6 in which said function is a weighted sum defined based on twoconstants W₁ and W₂.
 8. A method according to claim 7, comprisinglabeling possible values of voxel intensity by an integer index i andtheir respective frequencies by h(i), and writing the lower and upperintensities respectively as r_(low) and r_(high), wherein the weightedsum is given by:${{\theta_{RCLWV}\left( {W_{1},W_{2}} \right)} = {\max\limits_{r_{k}}\left( {{{\Pr\left( C_{1} \right)}{D\left( C_{1} \right)}W_{1}} + {{\Pr\left( C_{2} \right)}{D\left( C_{2} \right)}W_{2}}} \right)}},$where Pr(·) denotes we class probability$\left( {{\Pr\left( C_{1} \right)} = {{\sum\limits_{i = r_{low}}^{r_{k}}{{h(i)}\quad{and}\quad{\Pr\left( C_{2} \right)}}} = {\sum\limits_{i = {r_{k} + 1}}^{r_{high}}{h(i)}}}} \right),\lbrack,\rbrack$and D(C₁) and D(C₂) are given by:D(C ₁)=(μ₀ −μ _(T))², and D(C ₂)=(μ ₁ −μ _(T))², where${\mu_{T} = {\sum\limits_{i = r_{low}}^{r_{high}}{i \times {h(i)}}}},{\mu_{0} = {\sum\limits_{i = r_{low}}^{r_{k}}{i \times {h(i)}}}},{{{and}{\quad\quad}\mu_{1}} = {\sum\limits_{i = {r_{k} + 1}}^{r_{high}}{{i \times {h(i)}}.}}}$9. A method according to claim 1, wherein the steps (a) to (e) areperformed repeatedly, in each set of steps identifying at least onecorresponding landmark.
 10. A method according to claim 1, wherein steps(a) to (d) are performed to locate A, P, L and R landmarks, and whereinin step (a) the region of interest being defined within the AP plane;and in step (e) the most anterior and most posterior of the objectvoxels being taken respectively as the vertical coordinates of the A andP landmarks respectively, and the extreme horizontal components of theobject voxels are taken as the horizontal coordinates of the L and Rlandmarks respectively.
 11. A method according to claim 10 wherein step(c) comprises: performing at least one morphological opening operationon the binarized image obtained in step (b); and classifying one or morevoxels of the image(s) obtained by the opening operation(s) as objectvoxels or otherwise according to at least one criterion based ondistances in the image(s) obtained by the opening operation(s) andanatomical knowledge.
 12. A method according to claim 11 in which, priorto classifying the voxels, a maximum distance maxDSkull is obtained froma distance transform of the ROI.
 13. A method according to claim 11,wherein the classifying step comprises: restoring object voxels locatedfar from the skull, which were lost due to the morphological openingoperation(s); restoring object voxels around the boundaries lost due tothe morphological opening operation(s); and restoring object voxels lostdue to the partial volume effect.
 14. A method according to claim 1wherein steps (a) to (d) are performed to obtain an S landmark, themethod comprising: in step (a), defining the region of interest within avirtual plane obtained from Me a VPC coronal slice; and in step (e),identifying the position of the S landmark as the most superior point ofthe object voxels.
 15. A method according to claim 14 wherein step (c)comprises: performing at least one morphological opening operation onthe binarized image obtained in step (b); and classifying as objectvoxels one or more voxels of the image(s) obtained by the morphologicalopening operation(s) if they belong to eight voxels immediately adjacentto an object voxel and if their intensity value in the MR image ishigher than a value defined in relation to a second threshold.
 16. Amethod according to claim 1, wherein the set of steps (a) to (d) isperformed to identify an I landmark, comprising: in step (a), definingthe region of interest within a VAC plane; in step (e), defining the Ilandmark as the most inferior point of the object voxels.
 17. A methodaccording to claim 16 in which the threshold is obtained during apreceding process of locating an S landmark.
 18. A method according toclaim 16 in which, in step (c), (i) at least one morphological openingoperation, and/or (ii) at least one seeding operation, are performed onen the binarized image obtained in step (b).
 19. A method according toclaim 18 in which, in step (c), one or more voxels of the image(s)obtained by the morphological opening operation(s) which are notpresently classified as object voxels are re-classified as object voxelsif they are one of the eight immediate neighbours of an object voxel andif their intensity value in the MR image is higher than a value definedin relation to a second threshold.
 20. A method according to claim 16,wherein the left and right halves of the brain are treated separately,and the object voxels used to obtain the location of the I landmarkrelate to a selected half of the brain, the selected half of the brainhaving been selected based on a predefined criterion.
 21. A system forlocating one or more landmarks using an MR image of a brain, the systemincluding: an interface to receive data encoding the MR image; and aprocessor arranged to perform the following steps: (a) identifying aregion of interest with a plane of the MR image, the plane containingthe landmark(s); (b) binarising the plane of the MR image intoforeground and background voxels based on at least one thresholdselected using anatomical knowledge; (c) identifying a set of objectvoxels from the foreground voxels, the set of object voxels excludingvoxels which were only classified as foreground voxels due to proximityof cortical and non-cortical structures; (d) identifying object voxelsfrom the background voxels due to partial volume effect and/ormorphological erosion/opening; and (e) identifying the one or morelandmarks using the object voxels.